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Abstract 

Amplitude equations are derived that describe the spatiotemporal dynamics of cardiac alternans 
during periodic pacing of one- [B. Echebarria and A. Karma, Phys. Rev. Lett. 88, 208101 (2002)] 
and two-dimensional homogeneous tissue and one-dimensional anatomical reentry in a ring of 
homogeneous tissue. These equations provide a simple physical understanding of arrhythmogenic 
patterns of period-doubling oscillations of action potential duration with a spatially varying phase 
and amplitude, as well as explicit quantitative predictions that can be compared to ionic model 
simulations or experiments. The form of the equations is expected to be valid for a large class of 
ionic models but the coefficients are derived analytically only for a two-variable ionic model and 
calculated numerically for the original Noble model of Purkinje fiber action potential. In paced 
tissue, this theory explains the formation of "spatially discordant alternans" by a linear instability 
mechanism that produces a periodic pattern of out-of-phase domains of alternans. The wavelength 
of this pattern, equal to twice the spacing between nodes separating out-of-phase domains, is 
shown to depend on three fundamental lengthscales that are determined by the strength of cell- 
to-cell coupling and conduction velocity (CV) restitution. Moreover, the patterns of alternans can 
be either stationary, with fixed nodes, or traveling, with moving nodes and hence quasiperiodic 
oscillations of action potential duration, depending on the relative strength of the destabilizing 
effect of CV-restitution and the stabilizing effect of diffusive coupling. For the ring geometry, 
we recover the results of Courtemanche, Glass and Keener [M. Courtemanche, L. Glass, and J. 
P. Keener, Phys. Rev. Lett. 70, 2182 (1993)] with two important modifications due to cell-to- 
cell diffusive coupling. First, this coupling breaks the degeneracy of an infinite-dimensional Hopf 
bifurcation such that the most unstable mode of alternans corresponds to the longest quantized 
wavelength of the ring. Second, the Hopf frequency, which determines the velocity of the node along 
the ring, depends both on the steepness of CV-restitution and the strength of this coupling, with 
the net result that quasiperiodic behavior can arise with a constant conduction velocity. In both the 
paced geometries and the ring, the onset of alternans is different in tissue than for a paced isolated 
cell. The implications of these results for alternans dynamics during two-dimensional reentry are 
briefly discussed. 

PACS numbers: PACS numbers: 87.19.Hh, 05.45.-a, 05.45. Gg, 89.75.-k 
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I. INTRODUCTION 



Beat-to-beat changes in the T-wave portion of the electrocardiogram (ECG) that reflects 
the repolarization of the ventricles —a phenomenon known as T-wave alternans— jlL can 
be a precursor to life-threatening ventricular arrhythmias and sudden cardiac death [2]. T- 
wave alternans have been related to alternations of repolarization at the single cell level , 
thereby establishing a causal link between electrical alternans and the initiation of ventric- 
ular fibrillation. Repolarization alternans are characterized by a beat to beat oscillation of 
the action potential duration (APD) at sufficiently short pacing interval {4]. Experiments 
in both two-dimensional and linear strands of cardiac tissue, as well as ionic model 
simulations , have shown that the resulting sequence of long and short action po- 

tential durations can be either in phase along the tissue (concordant alternans), or can split 
into extended regions oscillating out of phase (discordant alternans). The latter case is of 
special importance since it can lead to conduction blocks 5] and the initiation of ventricular 



fibrillation [3|. Furthermore, alternans can provide a mechanism of wave breaks that sustain 
ventricular fibrillation j^, [^j . The dual role of discordant alternans in the initiation and, po- 
tentially, the maintenance of fibrillation provides the motivation to develop a fundamental 
understanding of this phenomenon. 



Experimental 



10] and theoretical studies Uj, ll2j have shown that spatially modulated 



patterns of alternans associated with quasiperiodic oscillations of APD can arise naturally 
in a ring of tissue. The terminology of spatially "discordant alternans" , however, originates 
from the observation of such patterns in periodically paced tissue jsj]. Their appearance 
in this context was first attributed to preexisting spatial heterogeneities in the electro- 



physiolo. 
present 



gical properties of the tissue 



HQ- 



Although heterogeneities may certainly be 



141 ]. and can affect the formation of discordant alternans |6j, numerical simulations 
of ionic models in homogeneous tissue jl], 0, Q, is| have shown that they are not necessary for 
their formation. Dynamical properties alone are able to produce the spatially heterogeneous 
distributions of APD observed in the experiments. 

Pioneering studies by Nolasco and Dahlen 16] and Guevara et al. [17] first explained the 
occurrence of alternans in a paced cell in terms of the restitution relation 

A p D n+l = J( DI n) ; (!) 

between the action potential duration at the n th + 1 stimulus, APD n+1 , and the so-called 



diastolic interval DI n . The latter is the time interval between the end of the previous (n th ) 
action potential and the time of the n th + 1 stimulus during which the tissue recovers its 
resting properties. The restitution relationship describes only approximately the actual beat 
to beat dynamics because of memory effects that have been modeled phenomenologically 



through higher dimensional maps I18|, [19|, 120|, 
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22 1. Moreover, recent theoretical studies 



have highlighted the role of intracellular calcium dynamics in the genesis of alternans [23|, [24 ] . 
In the present paper, we develop the amplitude-equation approach assuming that the simple 
relationship given by Eq. ([I]) holds. 

If the time interval between two consecutive stimuli is fixed, r = APD n + DI n = const. 
for all n, then Eq. ([T]) yields the map APD n = /(r — APD n_1 ) that has a period doubling 
instability if the slope /' of the restitution curve, defined by APD = /(DI), evaluated at 
its fixed point exceeds unity. As the slope of the restitution curve typically increases when 
the pacing period is decreased, alternans may develop beyond a critical value of the pacing 
rate. In extended tissue, the velocity of the activation wavefront also depends on DI, so 
the oscillations of APD induce oscillations in the local period of stimulation. This, in turn, 
acts as a control mechanism that tends to create nodes in the spatial distribution of APD 
oscillations, thus resulting in discordant alternans. 

In a previous paper [25(, we sketched the derivation of an amplitude equation that pro- 
vides a unified framework to understand the initiation, evolution, and, eventually, control 
of cardiac alternans 26l. 12?] . For one-dimensional paced tissue, the equation takes the form 



rd t a = aa — ga 3 — - / a{x)dx — wd x a + £ 2 <9^a, (2) 
A Jo 

where a ~ (APD n+1 - APD n )/2 measures the amplitude of period doubling oscillations 
in APD, so that nodes separating out-of-phase alternans region correspond to a = 0, r 
is the pacing period, a and g are coefficients that can be obtained from the map (TjQ), w 
and £ are lengthscales that depend on the strength of diffusive coupling between cells, and 
A is determined by the dependence of the conduction velocity (CV), denoted by c, on 
diastolic interval. The CV-restitution curve c(DI) is also known as the dispersion curve in 
the excitable media literature. It should be noted that time can be treated as a continuous 
variable in this framework because, even though the APD oscillates from beat to beat, the 
amplitude a of this oscillation defined above varies slowly over many beats close to the period 
doubling bifurcation. It is this slow spatiotemporal evolution that is described by Eq. (12]). 
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In this paper, we provide a detailed derivation of the amplitude-equation above and 
extend it to two-dimensional paced tissue. In addition, we extend this approach to treat the 
important case of one-dimensional reentry in a ring of tissue. For clarity of presentation, we 
will first derive the amplitude-equation for the ring geometry, and then show how to modify 
the boundary conditions on this equation to treat paced tissue in one and two dimensions. 
Although not treated in this paper, the amplitude equation formalism can be extended 
to higher dimensional maps (that include memory, or the effect of intracellular calcium), 
provided that the primary bifurcation is period doubling. 

The derivation of the amplitude-equation follows the general amplitude equation frame- 
work that has been widely used to study the evolution of weakly nonlinear patterns in 



non-equilibrium systems 



28] . In the present context of alternans, however, this derivation is 



made extremely difficult by the stiff nature of ionic models and the fact that the underlying 
stationary state (with no alternans) corresponds to a train of pulses. To surmount these 
difficulties, the derivation of the amplitude-equation proceeds in two steps. First, spatially 
coupled maps, which describe the beat-to-beat dynamics of APD oscillations, including the 
crucial effect of cell-to-cell coupling, are derived from the underlying ionic models. Second, 
a weakly nonlinear and multiscale analysis is used to derive amplitude equations from these 
spatially coupled maps. The weakly nonlinear nature of the expansion is valid close to the 
onset of the period doubling bifurcation. The multiscale nature of the expansion, which 
allows us to retain only the two lowest order spatial derivative terms in Eq. (j2j), is itself 
justified by the fact that the wavelength of spatial modulation of alternans is large compared 
to the length scale ~ £ over which the action potential of a cell is influenced by neighboring 
cells through diffusive coupling. 

The paper is organized as follows. In the next section, we introduce two ionic models and 
compute their action potential duration and conduction velocity restitution curves. These 
curves are used to relate quantitatively the ionic models and the amplitude equations. The 



models considered are the Noble model for Purkinje fibers [29] and a simple two-variable 
ionic model that captures basic dynamical properties of the cardiac action potential. This 
model has the advantage that it is simple enough to calculate analytically the restitution 
curves as well as all the coefficients of the amplitude-equation (j2j). 

Section II I II is devoted to the derivation of the amplitude equation in the ring geometry. 
This equation is then used to study the stability of pulses. The results are compared to 
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previous analyses based on coupled maps ll|, |l2j]. Courtemanche et al. UJ showed that 



the wavelength of spatial modulation of alternans is quantized in a ring geometry. In the 
same language that has been used recently to describe discordant alternans in paced tissue, 
this quantization condition corresponds, in the ring, to the existence of different modes of 
discordant alternans with an odd number of moving nodes separating out-of-phase regions of 
period doubling oscillations. The amplitude equation sheds light on the effect of cell-to-cell 
electrical coupling on the stability of these quantized modes of alternans. In particular, we 
find that this coupling lifts the degeneracy of an infinite-dimensional Hopf bifurcation so that 
the mode with one node is the most unstable, in qualitative agreement with the prediction 
of coupled maps that include a phenomenological description of cell-to-cell coupling 12]. 
Furthermore, we show that coupling can lead to quasiperiodicity even in the absence of CV- 
restitution. The effect of diffusive coupling on one-dimensional alternans dynamics in tissue 
has also been investigated in the context of a stability analysis of FitzHugh-Nagumo type 
models 30]. This analysis captures the fact that the onset of alternans in tissue can differ 
from a single cell. However, it does not provide a general framework to make predictions 
relevant for more realistic ionic models of cardiac excitation or experiments. 

We extend our theory to the paced case in section IIV1 The main result is that discordant 
alternans results from a pattern-forming linear instability that produces either standing, 
with fixed nodes, or traveling waves, with moving nodes that separate out-of-phase alternans 
regions. Whereas in the ring the wavelength of discordant alternans is determined by the 
ring perimeter, the wavelength in the paced case is independent of tissue size. In the latter, 
spatial patterns form due to the competing effect of CV-restitution, contained in the nonlocal 
term in Eq. (|2J), which tends to create steep spatial gradients of repolarization, and diffusive 
cell-to-cell coupling contained in the spatial gradient terms, which tends to smooth these 
gradients. As a result of this competition, the wavelength A of discordant alternans depends 
on the three fundamental length scales w, £, and A in Eq. (j2j), with different scalings 
A ~ (iuA) 1 / 2 and A ~ (£ 2 A) 1//3 for standing and traveling waves, respectively. We also 
extend the model to two-dimensional paced tissue. This extension is straightforward under 
the assumption that the propagation of the front is weakly perturbed by the oscillations in 
APD. Finally, we discuss briefly nonlinear effects and conduction blocks. 

Finally, in sections [V] and |VT] we present the discussion and conclusions. Appendix A 
contains additional details on the derivation of the restitution curves for the two-variable 
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ionic model. Appendix B, in turn, contains the details of the derivation of the integral kernel 
that couples spatially the iterative maps of the beat-to-beat dynamics. This derivation makes 
use of the Green's function for the diffusion equation to obtain an integral expression for 
the time course of the membrane voltage and hence the APD. 



II. IONIC MODELS AND RESTITUTION CURVES 

Our formulation is based on the action potential duration and conduction velocity resti- 
tution properties. To obtain the restitution curves, we simulate the ionic models in a one- 
dimensional strand of paced tissue. We consider the standard cable equation 

d t V = Dd 2 x V- (/ ion + I ext ) /C m , (3) 

with the membrane current Ji on , time in units of millisecond (ms), D = 2.5 x 10~ 4 cm 2 /ms, 
C m = 12/iF/cm 2 . The external current I ext models a sequence of stimuli applied at x = 
at a constant pacing interval r. We impose zero gradient boundary conditions on V at the 
two ends of the cable. 

In the rest of the paper, we will use two models for the ionic currents. The first is the 
original 1962 Noble model for Purkinje fibers [29j. The second is a simplified two- variable 
model with a triangular shaped action potential. This model is more tractable analytically 
for the purpose of deriving spatially coupled maps of the beat-to-beat dynamics and the 
length scales w and £ of the amplitude equation. For higher order ionic models like the Noble 
model, these coefficients need to be obtained numerically following a procedure detailed in 



Sec. IV. A. The two- variable model is defined by the total membrane current 25] 



^ = -[S + (l- S)V/V C ] + -h S. (4) 

This current is the sum of a slow, time-independent, outward current (similar to the potas- 
sium current in more realistic ionic models), which repolarizes the cell on a slow time scale 
r , and a fast inward current that is a simplified version of the sodium current. The latter 
depolarizes the cell in the fast time r a <C r . The inactivation of the fast inward current is 
regulated by the gate variable h, which evolves according to 

dh 1 — S — h 

dt = T_(l-S) + ST + , d 
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v c = 0.1 


e = 0.005 


r = 150 ms 


T a = 6 ms 


r_ = 60 ms 


t + = 12 ms 



TABLE I: Parameters used in the simulations of the two- variable model [Eqs. dU)-©] 

where r + and r_ control the time scales of inactivation and recovery from inactivation 
of this current. In this model the transmembrane voltage V is dimensionless, and S 
I + tanh((V — V c )/e)}/2 is a sigmoidal function. We consider the values of the parameters 
shown in Table [B Note that the product of the h and j gates in standard formulations of 
the sodium current is represented, in the present model, by a single gate h. Thus choosing 
r_ larger than r + produces the same effect of a j gate that controls a slower recovery from 
inactivation in comparison to inactivation controlled by the h gate. Fig. [1] illustrates the 
triangular shaped action potentials obtained with this model for stimuli spaced by 400 ms. 

We have simulated Eq. ([3]) using the forward Euler method, with a three-point finite 
difference representation of the one-dimensional Laplacian. We take a mesh size dx = 0.01 
cm, and a time step dt = 0.02 ms in the case of the two-variable model, and dx = 0.01 cm, 
dt = 0.05 ms, for the Noble model. The thresholds of the transmembrane voltage to define 
the APD and DI in the two- variable model and the Noble model are V = V c = 0.1 and 
V = —40 mV, respectively. The restitution and dispersion curves are computed by pacing 
Eq. ([3]) in a short cable by an S1-S2 protocol, i.e., the cable is paced at a large period for 
several beats (sufficient to achieve steady state APD), and then an extrastimulus (S2) is 
delivered at progressively shorter coupling intervals to vary DI (see Fig. [2]). The stimulus is 
typically applied over ten grid points. 

In the limit in which the sigmoidal function S = (1 + tanh((\^ — V c ) /e))/2 becomes a step 
function (e — > 0), it is possible to calculate both the APD- and CV-restitution curves of the 
two- variable model analytically (see Appendix |A} , provided that the APD is much larger 
than the time scale of inactivation of the fast inward sodium current, or APD ^> r+, and that 
D/c is much smaller than the wavelength of the action potential, or D/c < cAPD. Both 
conditions are satisfied for physiologically relevant conditions. In this limit, the restitution 
curve is found to be well approximated by the form 

APD ~ ^[1 - exp(-DI/r_)] - V c r . (6) 
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Accordingly, the onset of alternans is given by 

d APD r r_ 



dDI 



T„T- 



exp(-DI/r_) = 1, (7) 



or DI C ~ r_ log[(ror + )/(r a r_)], which yields DI C ~ 96.6 ms and APD C ~ 225 ms for the 
two- variable model parameters of Table [H 

In the same limit, the CV-restitution curve can be calculated considering a traveling pulse 
V(x,t) = V(x — ct), and matching the expressions obtained for V > V c and V < V c , in the 
wavefront and the waveback (appendix |Aj. This results in an implicit equation for c (cf. 
Eq. (IA20p ). that agrees well with the numerical results (see Fig. (2U). 

III. RING GEOMETRY 

The stability of pulses circulating in a ring of tissue was considered previously by Courte- 



manche et al. using generic restitution properties of the system as given by Eq. ([T]). This 
stability analysis was based on the idea to unravel the ring into an infinite line (then x G 1Z, 
instead of x G 1Z mod L). In this manner, the values of a given variable at a previous passage 
of the pulse through a given point x in the ring, and therefore, at a previous beat, can be 
identified with the value of that variable at the point x — L (i.e. APD n (x) = APD(x + nL)). 
This allows to drop the dependence of the different variables on the beat number n, and 
reformulate the maps as delay differential equations. Using this approach, Courtemanche et 
al. found that, when the slope of the restitution curve is greater than one, the pulses become 
unstable towards modes of propagation where the values of the APD and conduction veloc- 
ity vary in space around a mean value. Because of the periodic boundary conditions, the 
wavelength of these modes is quantized. Without CV-restitution, the ratio of the wavelength 
to twice the ring perimeter is l/(2n + 1) where n > in an integer. With CV-restitution at 
onset c'(DI c ) = dc/dDI(DI c ) ^ 0, this ratio becomes irrational and the nodes travel, giving 
rise to quasiperiodic motion at a given point. 



In the analysis of Uj all the modes were found to bifurcate simultaneously at the same 



ring perimeter, in an infinitely degenerate Hopf bifurcation. This degeneracy was shown to 
be broken [3] by the effect of diffusive coupling, which effectively modifies the restitution 
relation (OQ) and makes the mode with the longest wavelength (with a single node) bifurcate 
first. In [I2j, the effect of diffusive coupling was considered phenomenologically, assuming a 
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coupling of the APD in neighboring cells given by a simple gaussian kernel. Here, we derive 
the form of this kernel from an analysis of the two-variable model. The most interesting 
finding is that this kernel is generally asymmetrical due to the fact that parity (x — > —x) 
symmetry is broken by the propagation direction of the action potential; it only reduces 
to a symmetrical gaussian kernel in the limit of infinite conduction velocity. In effect, 
the way a given cell is influenced by its left and right neighbors is different because these 
cells are activated at different times by the propagating wavefront. The most interesting 
consequence of this asymmetry is that it can produce quasiperiodic motion even in the 
absence of CV-restitution, and modifies the frequency of quasiperiodic oscillations in the 
presence of CV-restitution. 

To study the stability of the pulses we will use an approach slightly different from that 
in We will convert the ring into a linear cable of length L with origin at x — 0, and 
index with n the number of times the pulse has traveled around the ring. Accordingly, the 
variables that characterize the beat-to-beat dynamics depend both on the beat number n 
and on the space variable < x < L; the diastolic interval is written as DI n (x) and the 
ring geometry imposes the periodic boundary condition DI™(x — L) = DI n_1 (x). Using this 
parametrization, the ring dynamics is governed by the equations 



where T n (x) is the time interval between two consequent activations at the same point. The 
first equation is just a kinematic equation stating that the period of stimulation at a given 
point is given by the time it takes a pulse to complete a revolution. The second equation 
reflects the restitution properties of the system, where we have included an asymmetrical 
kernel G{x' — x) that appears because of the diffusive coupling between neighboring cells. 
In general, the derivation of the kernel G{x' — x) is very difficult and we only derive it 
explicitly in Appendix [B] for our two- variable model. The calculation of the kernel involves 
the inversion of a nonlocal, implicit equation for the APD. Even though we cannot derive 
this kernel for a general higher order ionic model, this is not a serious limitation. Only the 
length scales w and £ in Eq. (T5]), but not the general form of the amplitude-equation, depend 
on a precise knowledge of this kernel. Moreover, w and £ can be calculated numerically for 
a general higher order ionic model using the procedure detailed in Sec. IV. A. 





(9) 
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A. Linear stability revisited 



Let us thus start reviewing the linear stability analysis of Eqs. jS]) and ©. Taking into 
account that T n (x) = APD n (x) + Dl n (x) one can write Eqs. (jSj), Q as a single equation 
for Dr +1 (x): 

rx A rp f poo 

DI " + w = L - L G(x ' ~ x)/[DI " MI dx '- (10) 

Then, a stationary solution satisfies DI* = r — /(DI*), where r = L/c(DI*) is the period of 
propagation of the pulse, and we choose the kernel G(x' — x) to be normalized, so 



/oo 
G(y)dy=l. (11) 
-oo 



We will also assume that the kernel has compact support and that it decays fast enough. 

Perturbing around the stationary solution in the form DI n (x) = DI* + a n e lkx 5D, then the 
stability of the steady state solution is dictated by the value of a. If \a\ > 1 an initially small 
perturbation will grow. Alternans occurs when a = — 1, and then the resulting instability 
gives a beat to beat change in the amplitude of APD. Inserting the former expansion in Eq. 
f TIU]) and linearizing, we obtain that the following equation must be satisfied 

a = ^(l-e~ lkL )-fG(k), (12) 

where we have defined the Fourier transform of the kernel 

G(k) = / G{y)^dy. (13) 



Here A = c 2 /(2c'), defined with c and d = dc/dDI evaluated at the bifurcation point, is a 
characteristic length scale associated with CV-restitution. 

Furthermore, from the condition DI n+1 (x — L) = DI n (x) we obtain the quantization 
condition e lkL = a, so Eq. (JI2"]1 can be rewritten as 

^-m^-f'^-m- (14) 

The onset of the instability is given by \a\ = 1 (implying that the bifurcating mode will 
have a real wavenumber k). This yields 

f' 2 \G(k)\ 2 + j^Lm[G{k)\ = 1. (15) 
11 



The quadratic equation can be solved to obtain the value of /': 

1 / Im[G(k)} 



f = [ ^ ± ^Im[G(kW/(Ak)^ + A\G(kWj . (16) 

The imaginary part of the kernel appears because of asymmetrical coupling. Due to the 
finite velocity of propagation of the pulses, different points are excited at different times, 
and the electrotonic coupling of a cell with its left- and right-neighbors will differ. We can 
assume that this is a small effect, and then: 

f ^ _J_ _ J HG(AQ] (17) 

\G(k)\ 2Ak\G(k)\ 2, 1 ; 

where we have taken the plus sign in Eq. ( fl6|) . corresponding to alternans (period doubling). 
The minus sign would correspond to a steady instability, not studied in this paper. 

Since the kernel decays on a spatial scale £ much shorter than the wavelength of the 
unstable modes of interest (~ L), the exponential in Eq. ffl3l) can be expanded in the form 
e iky r>j i _|_ fay _ (ky) 2 /2 ^ . It then follows that 

G(k) ~ 1 — iwk - £ 2 k 2 , (18) 
where we have defined the coefficients 

/oo 1 rco 

G(y)ydy, ? = - / G(y)y 2 dy. (19) 
-OO ^ J-OO 

For an arbitrarily complex ionic model, the form of the kernel G(y) cannot be calculated 
explicitly, and the coefficients w and £ 2 must be obtained from the numerical simulation 
of the ionic model as described in Sec. IV. A. For the two-variable model with a constant 
repolarization rate, however, they can be calculated explicitly (Appendix [B|, giving 

w = 2D/c, (20) 
f = (D x APD C ) 1/2 . (21) 

Eq. (1211) has the simple physical interpretation that the transmembrane potential V diffuses 
a length ~ £ in the time interval of one APD. Therefore, the repolarization of a given 
cell is influenced by other cells within a length ~ £ of cable. The imaginary part of G(k) 
appears because of the asymmetry induced by the propagation of the pulse. This asymmetry 
vanishes in the limit c — * oo, where all cells are activated simultaneously, consistent with Eq. 
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( 120]) . For typical values of the parameters in ventricular tissue (D ~ 10~ 4 — 10~ 3 cm 2 /ms, 
c ~ 10~ 2 — 10 _1 cm/ms, APD C ~ 100 — 200 ms), the length £ is of the order of millimeter 
(£ ~ 10 _1 cm), while w an order of magnitude smaller (w ~ 10~ 2 cm). In fact, w/£ ~ 
(/(cx APD C ) is the ratio between the diffusive coupling length £ and the wavelength of the 
pulse, which is typically an order of magnitude larger. 

Introducing expansion (1TB]) into Eq. (ITT]) , to lowest order in k£, we obtain: 

r = i-^+^ 2 . (22) 

Intercellular coupling shifts the onset of instability, which occurs for a value of the slope 
of restitution /' ^ 1. The slope at onset of alternans will be larger or smaller than one, 
depending on the relative importance of the different lengthscales of the system. The first 



12|. 



mode to bifurcate is the one with largest wavelength, as previously noted by Vinet 

In what follows, in addition to assuming that wk <C 1 and £k <C 1, we will also assume 
that the slope of the CV-restitution curve at the bifurcation point is small so that the length 
scale A is much larger that the typical wavelength of alternans modulation, or Ak ^> 1. All 
three conditions turn out to be reasonably well satisfied for the parameters of the Noble 
model. The last condition Ak ^> 1, however, is not fulfilled in the two- variable model 
because CV-restitution is not weak. In the next section, we will show that it is also possible 
to treat analytically the case where Ak is of order unity. 

Substituting Eq. (122]) into Eq. (TT4"j) . to lowest order we obtain 

ac^—l+iwk — — . (23) 
A A; v ' 

Both dispersion and asymmetrical coupling result in an imaginary part for a, that corre- 
sponds to quasiperiodic motion. 

If a = — 1, the instability is strictly period doubling. In this case, the condition e lkL = a 
implies k = (2n + l)ir/L, with n = 0, 1, 2, . . .. To obtain the wavenumber k corresponding 
to the value of a given in Eq. (123]) . we assume k = (2n + \)TxjL + q, where q is a small 
correction q 1/L. In this case the condition e lkL = a results into 

e ikL = ( _ 1)e ^ ^ + iqL + . . .) ^ _! + lw{ 2 n + X )Z[ _ lL +..., (24) 

L A(2n + l)7r 

Then, to first order in the correction, the wavenumber is 

k=(2n+l) i + mkvr,- w{ * n+v >T?- (25) 
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As k ~ L -1 , these expressions will be valid if w, £ <C L, and A 3> L. 

The main novel conclusion is that besides the correction to the wavenumber due to 
CV-restitution derived by Courtemanche et al. 11], there is another correction due to 
asymmetrical coupling. Therefore quasiperiodic motion can occur even with a constant CV. 
Furthermore, these two effects can balance each other if 

W = A(2ntl)% 2 ' (26) 
in which case the motion becomes strictly periodic even with a finite amount of CV- 
restitution. In general, these two effects will not exactly balance each other so that the 
motion will be quasiperiodic with a frequency determined by both the slope of the CV- 
restitution curve and the asymmetrical coupling strength ~ w. 



B. Derivation of the amplitude equations 

Starting from the maps (ED and (Q, it is possible to derive equations for the oscillations 
in period and action potential duration. To that end, we will consider the second iteration 
of the map (Q, and expand it for small values of the amplitude of oscillation. As the change 
in the value of the APD every two beats is small, the beat number can be treated as a 
continuous variable. Also, the dispersion relation for the conduction velocity c = c(DI) can 
be expanded for small oscillations of DI, obtaining from Eq. ([8]) the corresponding change 
in the local activation interval, which depends nonlocally on the oscillations of APD. The 
only non-trivial point in the expansion is the effect of the non-local electrotonic coupling in 
Eq. (Q. Since the kernel decays on a scale ~ £ shorter than the wavelength of modulation 
of alternans, it can expanded to obtain a local relation between the change in APD and its 
local first and second spatial derivatives. 

Close to the onset of oscillations we can write 

APD n (x) « APD C + 5 A + {-l) n a{x,t), (27) 
T n {x) « r c - St + (-l) n b(x,t), (28) 

where APD C and r c are the APD and the period of stimulation evaluated at the bifurcation 
point of the single-cell map (/' = 1), and we split the perturbations into steady and oscil- 
lating parts. Now the basic pacing period will be the traveling time of a pulse around the 
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ring, which in the absence of oscillations is given by r = L/c, and 5t = r c — r <C r c . Close 
to the bifurcation point, the steady correction to the APD is, to first order, 5 A = —5t/2. 
For the oscillating part, since the beat to beat oscillations are taken into account with the 
terms (— l) n , the amplitude of the deviations from the critical values, a(x, t) and b(x, t), vary 
slowly from beat to beat. We can therefore assume that a and b depend on a continuous 
time, defined through n = t/r. 

Let us first discuss what are the boundary conditions satisfied by a(x,t) and b(x,t). The 
transmembrane voltage obeys periodic boundary conditions V(L) = V(0), but, by definition, 
after a revolution of the pulse along the ring, the system goes into the next beat. Therefore, 
the values of the period and APD at x = L (that is, right before the end of the revolution) 
must equal those at x — 0, at the next beat (i.e. at the beginning of the next revolution). 
The same is true for the gradients. Then 

T n {L) = T n+1 (0), d x T n (L) = d x T n+1 {0), (29) 
APD n (L) = APD n+1 (0), d x APD n (L) = <9 x APD n+1 (0). (30) 

Using Eqs. ( 1271) . (1281) it is easy to see that, in terms of the oscillations in APD and period, 
the former boundary conditions become 

a(L) = -a(0), d x a(L) = -d x a{0), (31) 
6(L) = -6(0), d x b(L) = -d x b(0). (32) 

These boundary conditions imply that the pattern must have an odd number of nodes. 
Namely, letting a(x) ~ e* fcx , one obtains the quantization condition 

k = (2n+l)y, n = 0,1,2,..., (33) 

that obtained in ll| in the limit of zero CV-restitution slope. As we shall see, the corrections 
to the wavelength come from a phase shift in the slow scale associated with the quasiperiodic 
motion. 

The equation for the oscillations in period is easy to obtain. First, we can write Eq. (JHJ) 
in differential form 

dT n (x) = 1 1 = 1 1 r ^ 

dx " c[DI n (x)] c[m n {x-L)} c[DI n (x)] c[DI n - 1 (x)Y 1 ' 
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where we take into account that Dl n (x — L) = DI n_1 (x). Then, substituting expansions (J2"l 
and (1281) into the former expression, with DI n (x) = T n (x) — APD n (x), we obtain, at linear 
order 

^ = j[a(x)-b(x)). (35) 

To obtain an expression for b(x) we have to solve Eq. (1351) . subject to the boundary condition 
( 1321) . This results into: 

b(x) = 4- / e^'-^^a^')^' - — pTv- / e (a; '- x)/A a(a; , )^ / . (36) 

A Jo A(l + e- L / A ) jo 

In order to be able to get an analytical expression for the shift in wavelength, we will take 
the limit of small dispersion (L/A <C 1), that will also allow us to compare with the results 



m In this limit, to first order, the exponentials in Eq. (1361) can be neglected (equivalent 
to neglecting the term proportional to b(x) in the right hand side of Eq. (135]) ). Then, Eq. 
( 1361) becomes 

b( x ) = \ ! X a{x')dx' -—I a{x')dx'. (37) 
A jo 2A Jo 

The above equation works well for the Noble model. For the two-variable model, however, 
it is not strictly valid since L c /A > 1. In this case, as we discuss below, the full equation 
(|36|) should be used to obtain a good agreement with the simulations of the cable equation 
(J3]). For clarity of exposition, we use Eq. ( !37l) in what follows and state the result later for 
the full equation ( 1361) in the context of a quantitative comparison of the two- variable model 
and cable equation simulations. 

Next, in order to derive an evolution equation for the amplitude a(x,t), we notice that, 
after two consecutive beats, the APD becomes 

ApD „ +2 = ApDc + (_!)n+2 a ^ t + 2r ). (38) 

Assuming that the APD varies slowly every two beats (which is the case close to the period 
doubling bifurcation), we can expand a(x,t + 2t) ~ a(x,t) + 2rda/dt, so 

ApD n+2 = ApD n + (_ 1 )n 2r ^_ > (39) 

But, expanding Eq. ([9]) 

/oo 
G{y)f[m n (y + x)} dy ~ J[DP(x)] - wf'd x Dl n (x) + £ 2 f ^Dr(x), (40) 
-oo 
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we can also write APD™ +2 as 

ApD n+ 2 = /[ T «+i_/( r «_APD™)+wf9 ;E Dr-e 2 f^Dr] 

-wfd x (T n+1 - APD n+1 ) + £ 2 fd 2 x {T n+l - APD n+1 ), (41) 

where we have used the definitions in Eq. (1T91) for w and £ 2 . Then, using Eqs. ( l2~Tj) - (l28l) . 
expanding to first order in w and £ 2 , and taking into account that f = 1, since we are 
expanding around the period doubling bifurcation, we have 

ApD n+2 = jpn+l _ f( T n _ A p D ")] + ^DI" - £ 2 a£DF 

-wd x [T n+1 - f(T n - APD n )] + £ 2 <9 2 [T n+1 - f(T n - APD n )] 
= f[ T n+i _ j( T n _ APD n)] + (_!)« \_ w ^d x a - 3d x b) + £ 2 (2<9 2 a - 3d 2 b)] (.42) 

In the following we will neglect the term d x b compared to d x a, since from (1371) we have that 
b ~ a/(Ak) a, in the limit of small dispersion (i.e. weak CV-restitution) . Strictly, the 
amplitude equation we derive is asymptotically valid close to the bifurcation point if the 
following scaling relations are satisfied: a/APD c ~ e 1 / 2 , b/a ~ l/(A/c) ~ e, £k ~ e 1 ' 2 , and 
wk ~ e, where e ~ St/t c is a dimensionless measure of the distance from the bifurcation 
point. 

Equating now Eqs. fl39l) and fT42|) . and expanding the second iteration of the map, we ob- 
tain the final expression for the evolution of the oscillations of APD, in the limit considered, 

rd t a = aa — wd x a + £ 2 <9 2 a — ga 3 — b, (43) 

where a = f"(r — r c )/2 is of order e, g = /" 2 /4 — f'"/6, and all the derivatives are evaluated 
at the bifurcation point. These coefficients can be calculated from the curves in Fig. [21 
We find the bifurcation in the Noble model to be slightly subcritical, so Eq. (H3l) must be 
expanded to fifth order in this case 31(. For simplicity of exposition, in the following we 
will focus on the supercritical case. When dispersion is not small, one should keep higher 
order terms in Eq. (143]) involving the oscillation in period b. In that case, direct simulation 
of the original coupled maps (JHJ)- (ED is probably more appropriate, if the goal is to achieve 
good quantitative agreement with ionic model simulations. 

Substituting the expression for b(x) into Eq. ( 1431 . we obtain the final amplitude equation 
in the ring geometry 

1 r x 1 fL 

rd t a = aa — ga 3 — wd x a + £ 2 <9 2 a — - / a(x')dx' + — / a(x')dx', (44) 

A jo 2A jo 
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that must be solved with the boundary conditions a(L) = — a(0), d x a(L) = —d x a(0). 

Now, we can consider again the linear stability problem, within our amplitude equation 
framework. To do that, we write a(x) ~ e ikx+sit^ w j^h k given by (I3"3"j) . and Q complex 
(Q = Q r + iQi). Separating into real and imaginary parts 

f 2 

T^ r = a-i 2 k 2 = a- ^r(2n+ 1) 2 tt 2 , (45) 

L 2 

rQi = -L - wk= --, — ■ tt ^(2ra + l)7r. (46) 

Afc A(2n+l)f L l ; v ; 

The onset of instability occurs when Q r = 0, which results in a = £, 2 k 2 , equivalent to 
condition fl22|) when w/A <C 1. Again, we see that intercellular coupling lifts the degeneracy 
in the onset of the different modes, since the growth rate depends on the wavenumber of 
the mode. Clearly, the fastest growing mode is that with the largest wavelength, which 
corresponds to the mode with a single node. Intercellular coupling also affects the frequency 
of the quasiperiodic oscillations. 

To compare Eqs. (J45]), f|46|) with the results from the linear stability of the maps in the 
previous section, we can factor out the rapid oscillations, so a n = (— l) n f3 n = (—/?)", and 
(3 n = e Qt = e TQn , by the definition of the slow time t = nr. Then, at onset of the instability, 
Q = iQi, and 



a = -(3 = -e irn > ~ -1 +iwk- 4r, (47) 

Ak 



which is exactly the expression we had obtained before (cf. Eq. (|23|) ). 

We have checked these results with simulations of the two-variable model. As expected, 
the first mode to bifurcate presents only one node, and the oscillations are quasiperiodic (see 
Fig. [3]). The onset of alternans appears in the simulations for L c = 5.11 cm, and the phase 
speed of the bifurcating mode is v — —Qi/k = —1.81 • 10~ 3 cm/ms, which can be computed 
from the position of the node as a function of time (see Fig. [3]). Using the restitution curves, 
one obtains that /' = 1 when r c = 321.5 ms, for a conduction velocity of c = 0.0161 cm/ms. 
This yields the prediction L c = 5.17 cm. This prediction, however, neglects the stabilizing 
effect of the diffusive coupling. Setting k = n/L, and using f|45|) with the parameters of the 
two- variable model, the predicted critical ring length that includes this effect is L c = 5.13 
cm, which is in good agreement with the value L c = 5.11 cm in the numerical simulations 
(Fig. H|). Then, setting L = L c , the critical wavenumber is k c = tt/L c = 0.608 cm -1 , and 
Eq. (l4"6]) gives fli = 1.38 • 10~ 3 ms -1 at the onset of instability, which results in a velocity of 
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the nodes (equal to the phase velocity) v = —Qi/k = —2.27 • 10~ 3 cm/ms. The discrepancy 
with the simulation value is due to l/(A/c c ) corrections. Starting from the full equation fl35l) . 
which does not assume Ak c ^ 1, it is possible to obtain the modified expression for the 
frequency 

Ak 

Ttti = T+JaW ~ wt (48) 

This expression gives v = — 1.85-10 -3 cm/ms, in almost perfect agreement with the numerical 
results. 

Close to onset, the bifurcating solution is therefore generally in the form of a traveling 

wave 

a(x) = ^(Be i{kx+n ^ + c.c), k = tt/L. (49) 

Substitution this form in the full nonlinear amplitude equation for the ring ( ]44"1) and balanc- 
ing separately real and imaginary parts, we obtain at once that the bifurcation is supercritical 
with a traveling wave amplitude 



|B| = ii- — r- 1 - (50) 

Given the expected universal validity of the amplitude equation (in the limit of small dis- 
persion A — > oo), this result implies that the bifurcation to alternans in the ring will always 
be supercritical if the quasi zero- dimensional bifurcation to alternans in a periodically stim- 
ulated small tissue patch is supercritical, which occurs when g > as in the two-variable 
model. (Note that it is important to distinguish a small tissue patch from an isolated cell 
that has a different restitution curve due to the absence of diffusive coupling.) We have 
checked numerically that the bifurcation is supercritical in the two- variable model (Fig. HI). 
The amplitude predicted by Eq. (50), however, is about 60% higher than the numerically 
computed amplitude shown in Fig. HI This discrepancy can be attributed to l/(Ak c ) cor- 
rections that are not small in the two-variable model due to steep CV-restitution. It should 
be possible to improve the agreement by including these corrections in the weakly nonlinear 
analysis. Conversely, if g < 0, the bifurcation in the ring is subcritical and the saturation 
amplitude is generally determined by higher order nonlinear terms in the amplitude equa- 
tion, which can be computed analytically if the bifurcation is only weakly subcritical as in 



the case of the Noble model 



3l|. 
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IV. PACED TISSUE 



The main difference between the ring and the paced case comes from the role of the 
boundary conditions. While in the ring the periodic boundary conditions result in a quanti- 
zation condition for the unstable modes, such a condition is absent in the paced case. Thus, 
the selected wavelength in this latter case must be related to some intrinsic length scale 
of the system. We show in Fig. [5] simulations of both the Noble and two-variable mod- 
els. We have done simulations in long cables to highlight the striking spatial regularity of 
out-of-phase domains of alternans. These patterns suggest that the formation of discordant 
alternans is caused by a finite-wavenumber linear instability of the basic underlying state, 
similar to those encountered in other physico-chemical systems, such as Rayleigh-Benard 
convection, Taylor-Couette flow, etc 28]. The APD oscillations obtained with the Noble 
model resemble a standing wave (Fig. [5b), and those for the two- variable model a traveling 
one (Fig. (5bl). Thus, at a given point in the tissue, the oscillations in APD are periodic in 
the first case, and quasiperiodic in the second. We will see how, within the formalism of the 
amplitude equations, these length scales appear naturally. 

To derive the amplitude equations for paced tissue, we must write the equivalent of Eqs. 
(JHJJ-Q for this case. Now, the period at a given point will be given by: 

T " (X) =T + i c[Df(x')] ~ I c[Dr- 1 (x')Y (51) 
where r = T n (0) is the period of stimulation at the pacing point. This simply means that 

the period of stimulation at a given point in the tissue is equal to the period of stimulation at 

the pacing point plus the difference in arrival time between two consecutive pulses. Notice, 

however, that in differential form this equation is the same as in the ring (cf. Eq. (1341) ) . 

To complete the system, we can use again Eq. but with a note of care. In fact, in the 

derivation in the ring, we use translational symmetry to expand the kernel. In the paced 

case this symmetry is broken and one should, in principle, calculate the kernel for the finite 

system. When the decay rate of the kernel is fast this does not seem to be necessary. Eqs. 

(I5ip -(I9]L supplemented with non-flux boundary conditions for the oscillations of APD, give a 

remarkably good agreement with simulations of the ionic models, for all the lengths of tissue 

considered in this paper. It is interesting to notice that the APD itself does not satisfy these 

boundary conditions, as rapid spatial variations in APD ("blips") appear at the two ends 

of the cable due to the non-flux boundary conditions for V. The amplitude a of alternans 
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obtained by taking the difference of APD between two consecutive beats, however, does 
satisfy very well the non-flux boundary conditions. This can be checked from numerical 
simulations of the ionic models (see Ref. 26|). 

It follows that the amplitude equations in the bulk remain the same [cf. Eq. (|43p ], and 
we only have to modify the boundary conditions. The condition that the period must be 
equal to the pacing interval at x = 0, T n (0) = r, translates into 6(0) = 0. Then, Eq. (|35|) 
can be solved for b(x) to give 

b(x) = j£ e {x '- x)lk a(x')dx'. (52) 

As in the previous section, we will assume that the CV-restitution curve is shallow at the 
bifurcation point, so A is much larger than the wavelength of the pattern (A ^> 2n/k). If this 
is the case then the exponential in (152"j) can be neglected and the former relation becomes 

b{x) ~ - f a(x')dx'. (53) 
A Jo 

Introducing this into Eq. (l4"3j) . we arrive at the final expression 

1 f x 

rd t a = aa — ga 3 — — J adx' — w d x a + £ 2 <9^a. (54) 

With the help of equation (154"]) the formation of the nodes is easy to explain. Let us first 
neglect the influence of the electrical coupling between cells on the APD (that later on will 
be shown to be crucial), and consider the equation 

rdta = aa — ga 3 — b. (55) 

Starting at x = with a constant amplitude of oscillation for the APD do, the oscillation in 
the period b will increase along the tissue [cf. Eq. ( l53l) ]. effectively decreasing the value of a, 
until it crosses zero and goes into the branch with opposite phase (see Fig. [6]). Thus, the sys- 
tem tends to create discordant alternans spontaneously, through CV-restitution. However, 
without the derivative terms the gradients become ever steeper, tending to a discontinuous 
limit as t — ► oo. To see this, we can differentiate Eq. (153|) . to obtain the steady state equa- 
tion for the slope of the oscillations in APD d x a = a/A(a — 3ga 2 ), that can be integrated to 



obtain x = A[crln(a/a ) — 3g(a 2 — af ) )/2], with ao = a(0) = \Jo-/g. When a 2 (x ) = cr/Sg the 
slope becomes infinite, denoting the formation of a singularity at xq. The singular behavior 
originates from the use of the local APD-restitution relation [Eq. ([T|)] that allows two nearby 
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cells to have different APDs, whereas this is prevented in the cable equation by the spa- 
tial diffusive coupling. The distance between singularities can be obtained integrating Eq. 
(I55|) from a 2 = 4er/3<? (the value after the singularity) to a 2 = cr/3g (the value at the next 
singularity). Then, we obtain the wavelength (as twice the distance between singularities): 
A = 2Acr[3/2 — ln(2)], that does not correspond to the right length scale for this problem 
(compare Figs. [6] and [5]). Once the derivative terms are included, the agreement between 
the ionic models and the amplitude-equations becomes very good. To check this point, we 
have derived the coefficients w and £ 2 from the restitution curves of the Noble model (in the 
two- variable model their values are given by Eqs. ( [20]) and ( T2~B ). 

A. Measurement of the coefficients w and £ in the amplitude equation. 

To measure the coefficients w and £ we will consider a case where tissue is paced at a 
constant period, for values close to, but beyond, the onset of alternans. If we plot the map 
APD n+1 in terms of DI n corresponding to those simulations (Fig. [7J we obtain a restitution 
curve that differs from the one obtained using a S1-S2 protocol, where DI is spatially uniform. 
During discordant alternans the electrotonic currents generated by the spatial modulation 
of DI along the cable modify the cable restitution curve computed in Appendix |A] for a 
spatially uniform DI, which is equivalent to the one obtained using a S1-S2 protocol, fsiS2 
(see Fig. [2]). The cable restitution curve with a spatially modulated DI, which we define as 
fmod, can be different on even and odd beats during discordant alternans because the sign 
of the spatial gradient of DI in the nodal region alternates from beat to beat. If the system 
presents no memory, then the difference between f mo d and fsiS2 is due to the presence of 
gradients of DI. Expanding the kernel in Eq. (Q, we can write: 

APD" +1 (x) = / slS2 [Dr(i)] - wd x DI n (x) + e 2 ^DI"(x) = f mod [m n {x)}. (56) 

Let us now consider a point in space x* where there is a node in DI, so DI n+1 (x*) = DI n (x*) = 
DI*. Then, the S1-S2 restitution curve gives the same value at two consecutive beats, so the 
split in the restitution curve in Fig. [7] will be due to the gradients. Given that the alternans 
profile in Fig. [T] is, to a good approximation, sinusoidal, then at that point it is also satisfied 
that c^DI n (x*) ~ d|DI n+1 (x*) ~ 0. Subtracting the APD at consecutive beats we have: 

AAPDi = f mod [DI n+1 (x*)) - / morf [DF^*)] = w[d x Dl n (x*) - ^Df +1 (x*)]. (57) 
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Model 


A 


w 




\heorl 4 


-^sim/4 




Noble 


49.1 


0.045 


0.18 


2.33 


2.6 


2.75 


Two- variable 


3.55 


0.031 


0.235 


1.33 


1.1 


1.15 



TABLE II: Values of the parameters in the amplitude equation ([2]) associated with different length 
scales, and comparison of the simulated and theoretical wavelengths for the unstable modes in the 
paced case (all distances in cm). The coefficients w and £ are obtained analytically from Eqs. (|20|) - 
(|2ip in the two-variable model, and numerically from the APD-restitution curves as explained in 
the text, for Noble model. The values of A are obtained from the CV-restitution curves computed 
numerically, in both cases. The wavelength Xtheor is predicted from Eq. (f67l) . for Noble, or Eq. 
(|65p . for the two- variable model (see text for details), while \ s i m is obtained from simulations of 
Eq. ([3]) in a long cable. These wavelengths give also a good estimate of the minimal length of 
tissue at which a node appears. 



Then, we have 



w 



AAPDi 



d x Bl n (x*) - d x Dl 



n+l. 



X 



(58) 



To calculate the coefficient £, we just consider the point of tissue x max where f m od has a 
local maximum, corresponding to an antinode of alternans. At that point d x Dl n (x max ) = 0, 
and then: 

f mod [Dl n (x max )} = fsis2[VI n (x max )} + ed 2 x DI n (x max ). (59) 

In this case £ can be obtained from the difference between the S1S2 restitution curve 
and the one obtained during discordant alternans. Defining AAPD2 = fm d[Dl n (x max )] — 
fsiS2{Ur(x max )}, we have: 



AAPD, 



(60) 



d 2 x m n {x max )' 

In Table ILT1 we show the values of the coefficients calculated in this manner. The comparison 
between simulations of the amplitude equation (!54j) using these coefficients and the Noble 
equations is very good, as shown in Fig. [HJ 



B. Linear stability analysis 

Although nonlinear effects determine the saturated amplitude of alternans, the spacing 
between nodes of discordant alternans, and the velocity of the nodes in the case of traveling 
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waves, are well predicted by linear stability theory close to the bifurcation. 

The genesis of discordant alternans can be understood by computing the linear stability 
spectrum of the spatially homogeneous state (a = 0). We have calculated this spectrum 
numerically for different values of L, and compared it with the analytical results obtained in 
the large L limit. The main result is that the wave pattern can emerge from the amplification 
of either a unique finite wavelength mode, which yields a stationary pattern, or from a 
discrete set of complex modes that approach a continuum in the limit of large L, and yields 
a traveling pattern (see Fig. [5]). There is indeed experimental evidence for both stationary 
jsl and traveling Q| waves. 

In the large L limit it is possible to obtain analytical expressions for the wavenumber 
and onset of the waves, using the dispersion related associated to Eq. (|54|) . Differentiating 
this equation with respect to x and letting a(x,t) ~ e ^x+nt/r^ DO th Q = f2 r + iVti and 
k = k r + iki complex, yields at once 

Q — a — i 2 k 2 - i 



wk — 7-7- 

Ak_ 



(61) 



Except for reflection symmetric systems (i.e. system invariant under the transformation 
x — > — x), or with periodic boundary conditions that impose k G Tie, the wavenumber will 
be complex. An instability occurs if there exists a complex wavenumber k for which Q r > 0. 
It may happen, however, that the unstable mode has a non-zero group velocity dfl/dk ^ 0. 
Then, the instability will grow in the frame moving with the group velocity, but decay at a 



fixed position. This is the signature of a convective instability [28| where perturbations are 
transported as they grow, similarly to, for instance, Taylor-Couette vortices developing in 
an axial flow j^. In such a situation, patterns are only transient in a finite system (unless a 
constant forcing is provided), since they disappear through the boundaries in finite time. A 
sustained instability develops if the group velocity of the unstable mode is zero dfl/dk = 0, 
in which case the instability grows at a fixed location in space. This is known as an absolute 
instability (Fig. [5]). The condition 

«.o— aC*-*.-^ (62) 

yields a prediction for the complex wavenumber corresponding to the onset of absolute 
instability. In the limit w — > 0, it becomes simply 

k = ± yl l __ f63 ) 

2(2eA) x / 3 2(2eA) 1 / 3 ' 1 ' 
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and corresponds to modes growing exponentially in space. This mode is evident in Fig. [5H, 
where the amplitude grows away from the pacing point before saturating due to nonlinear 
effects. There exists a third solution of Eq. (I62p that gives a purely imaginary wavenumber 
(an exponentially decaying mode without sinusoidal part). It can be checked, however, that 
this mode belongs to the spectrum associated with the equation resulting from differentiating 
Eq. (1541) with respect to x, but not to Eq. (1541) itself. 

Substituting expression (1631 into Eq. (ISTl) . we deduce that the threshold of absolute 
instability occurs when 

a^ = (3/2)(£/2A) 2 / 3 , (64) 

with a pattern of wavelength 

A = ^|(2£ 2 A) 1/3 . (65) 

This wavelength agrees well with that observed in simulations of the two-variable model 
(see Table HTjl . The frequency is, in this approximation, fi; = (3>/3/2)(£/2A) 2 / 3 , resulting 
in Qi ~ 0.27 for the two-variable model, and the pattern travels with phase velocity v p h = 

We have confirmed the validity of the former analysis by numerically solving the linear 
eigenvalue problem associated with (}54"1) . To obtain the linear spectrum for a given value of L, 
we have linearized and discretized in space Eq. (1541) . using a finite difference representation 
of the derivatives and the trapezoidal rule for the integral. Looking for exponentially growing 
or decaying solutions di(t) ~ aie nt ^ T , we obtain a set of N coupled linear algebraic equations 

w £ 2 dx 1 

Qcii = aa i -^^(ai +1 -a i . 1 ) + — (ai +1 +ai- 1 -2ai)-—^2-(a j +a j+1 ), i = 1...N, (66) 

where L = Ndx, and we typically use dx = 0.05 cm. The non-flux boundary conditions 
now become a = a 2 , and a^+i = un-i- The resulting eigenvalue problem is then solved 
for the complex growth rate Q = Q r + iQi, and the corresponding eigenmodes. These will 
be either real (with fij = 0), in which case we will talk about stationary modes, or come 
in complex conjugate pairs, resulting in traveling waves. In order to compare with the 
analytical predictions, we calculate the wavenumber from the number of nodes n of the 
eigenmodes (of both their real and imaginary parts), as k r ~ 7m/ L. 

From Fig. [10] it can be seen that the spectrum consists of a continuous branch, plus an 
isolated eigenvalue, whose origin is probably due to the boundary conditions. This can be 
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motivated by noting that cos k r x is an exact eigenvector of Eq. (15^1) linearized around a = 
that satisfies d x a = at the two cable ends. Substituting that solution into Eq. fl54|) . we 
obtain f2j = 0, Q r = a — £ 2 /c 2 , and a wavelength A = 2n/k r given by 

A = 2n(wk) 1 ' 2 . (67) 

This is in good agreement with the wavelength of standing waves observed in simulations of 
the cable-Noble equation (Table HT1). Although this exact mode only exists for given specific 
values of the length of the system (when L is an integer multiple of A/2), the computed 
spectrum shows that it also persists for arbitrary values of L. The onset of standing waves 
then is predicted to occur at 

a th = ?/(wA). (68) 

The continuous spectrum ends in the branch points given by the absolute instability, 
with both k and fl complex, corresponding to traveling waves. Then, the bifurcating modes 
can be either stationary or traveling waves, depending if it is the isolated eigenvalues, or 
the branch points, that bifurcate first. In the limit w — > the growth rate of the standing 
wave fl r — > — oo, and the resulting pattern will always be that of traveling waves. There is, 
therefore, a critical value of w below which the isolated eigenvalue presents a lower growth 
rate than the branch points, and complex modes that grow exponentially at large x are 
the most unstable. This is the case for the two- variable model (see Fig. fTOT) . where the 
modes pass from being convectively to absolutely unstable, as the pacing rate is decreased 
(Fig. [TT|) . From Eqs. (}6"lj) and fl68l) it can be deduced that traveling waves are favored over 
stationary waves when A = c 2 / (2c') <C £ 4 /u> 3 , and hence for strong CV-restitution, and vice 
versa for weak CV-restitution. 

Thus, our results demonstrate that the formation of discordant alternans is crucially 
affected by the effect of electrical coupling (diffusion) on repolarization, in addition to resti- 
tution and dispersion. Dispersion is responsible for the formation of nodes and spatial 
gradients of APD that steepen with time. Diffusion, in turn, tends to spread the APD 
spatially, and also induces a drift of the pattern away from the pacing site that is induced 
by the more subtle gradient term (— wd x a) in the amplitude equation. When dispersion 
is sufficiently weak, it may be balanced by drift and produce a stationary pattern. In the 
opposite limit, the tendency for dispersion to form steep gradients of APD is balanced by 
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the spreading effect of diffusion. Nodes then travel, cyclically disappearing (appearing) at 
the pacing (opposite) end of the cable. 

The length scale for the standing waves in the paced case is the same as that for the 
strictly periodic motion in the ring [see Eq. (|26p ]. In the former case, however, the non-flux 
boundary conditions pinned the node (imposed Qi = 0), and that was the only possible 
wavelength for a bounded state. The only possibility of having a different wavelength, and 
a traveling node, was to pick up an exponential part (k complex), which is precluded in 
the ring by the boundary conditions. There, the wavelength is selected by the boundary 
conditions, and depends on the length of the system. This, in turn, selects the velocity 
of the node through Eq. (146 p . that in general will be different from zero. Comparing the 
critical length scales in both cases, we see that for the mode with largest wavelength in 
the ring (n = 0), this results in a critical tissue length of L c = A/2 = 5.78 cm, for Noble, 
and L c = A/2 = 5.17 cm for the two- variable model. Thus, the length scale in the ring is 
comparable with the one in the paced case for the Noble model, while for the two-variable 
model it is significantly larger (see Table Hill . 

C. Stability diagrams 

To compare with previous numerical and experimental results we have calculated the 
different solutions of the Noble and two- variable models as a function of cable length L and 
pacing period r. As can be seen in Fig. [12] these results are in good agreement over a wide 
range of L and r with those obtained simulating the amplitude equations. In particular, 
discordant alternans appear as soon as the cable is long enough to accommodate roughly 
a fourth of a wavelength, as given by expressions flBTj) and (165]) (see Table fld])- Several 
comments are in order regarding our results. 

(i) The transition from concordant to discordant alternans occurs as the length of the 
tissue is increased, but the transition line remains nearly constant with increasing pacing 
rate. This seems to contradict the results obtained in [3], where a transition from concordant 
to discordant alternans was reported as the pacing rate was increased. In typical models, the 
slope of CV-restitution is small for large values of DI, but increases very rapidly for small 
values of DI (see, for instance, the CV-restitution curve of the Noble model in Fig. [2b). 
Thus, as the pacing rate is increased, the amplitude of alternans grows thereby engaging 
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the steep part of the CV-restitution curve and causing a transition to discordant alternans. 
This does not occur in the Noble model since there is a reverse period-doubling bifurcation 
at larger pacing rates by which the steady state regains stability, and the oscillations never 
get to grow enough, but we have checked this point with simulations of the Beeler-Reuter 
equations (not shown here). Within our amplitude equations this effect can not be captured 
since we are expanding around the period-doubling bifurcation point, where the slope of the 
CV-restitution is fixed. However, it can be observed with the coupled map equations. 

(ii) Discordant alternans appear directly from the rest state. This was already observed 
in simulations of the Beeler-Reuter model p|, for lengths of tissue L > 5 cm. As can be 
seen in Fig. 5 of that paper, the number of nodes increases with the size of the tissue. This 
is also the case in our simulations of the Noble model (Fig. \5§ where, close to onset, the 
oscillations of APD adopt a sinusoidal form, with a well defined wavelength. A similar result 
was observed in |6j], where the Luo-Rudy model was modified to obtain several APD and 
CV-restitution curves. In agreement with our results, when the slope of CV-restitution was 
large at onset, a direct transition to discordant alternans was observed, while, in the case of 
a shallower curve, the primary transition was to concordant alternans, discordant alternans 
only appearing when the oscillations in APD grew enough to explore the steep part of the 
CV-restitution curve. 

(iii) The onset of alternans in an extended tissue is delayed with respect to that of the 
single cell (as given by Eq. (CQ)). This restabilizing effect is due to a combination of diffusive 
coupling effects and CV-restitution (see Eqs. (|68l)-(l64l)). that induces oscillations in the 
period of stimulation that effectively act as a control mechanism. In Eq. (15^1) this effect is 
produced by the integral term. As is readily seen in Fig. [12] from the diagrams of the Noble 
and two-variable models, the larger the slope of CV-restitution is at the critical point /' = 1 
(the smaller A is), the larger the value of the restabilization. 

(iv) The nodes separating regions oscillating out of phase may be stationary (Noble 
model), or traveling (two- variable model). Since this distinction occurs already at onset, 
and there is no transition from one case into the other as the pacing rate is increased, the 
selected type must depend on the specific ionic model considered. There is experimental 
evidence for both types of patterns [3, |5|]. In real tissue, however, other effects not considered 
here, such as gradients of restitution properties, could be important, and affect the motion 
of the node. In the experiments in [5] the tissue was paced from both ends, with identical 
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results, which seems to suggest that gradients are not important in that case. In other 
situations, however, the gradient terms could become important, and create a drift that 
opposes or reinforces the effect of dispersion. 



D. Nonlinear state: conduction blocks and front solutions 

Nonlinearities are mainly responsible for saturating the amplitude of the linear modes, 
although both the wavelength and evolution of the pattern generally vary with distance 
from onset. Furthermore, at high amplitudes conduction blocks can be produced, which 
are the main mechanism for the induction of reentry. Recently, systems of coupled maps 



simi 



33, 



ar to 



qs. ([T]) and © have been used to study the appearance of conduction blocks 
351 ] . They are often produced right after the first node in the APD oscillations 
33j, thus stressing the importance of discordant alternans for the induction of reentry. 
Furthermore, the_position of the wave block has been shown to vary if the system presents 
traveling nodes [34j]. 

Although strictly valid only close to onset, we will use our amplitude equations to get 
an idea on where conduction blocks can be formed. A conduction block occurs wherever 
the APD is large enough, so DI = r — APD < DI m j n , being Dl min the minimum diastolic 
interval necessary for propagation. Let us first neglect the diffusive coupling. Then, using 
Eq. (155]) . which is the amplitude equation corresponding to the maps flTJ and flS]), it is easy 
to show that the maximum value of the APD occurs just after the discontinuity (see Fig. 
[6]). Considering that the system has reached steady state, and taking the spatial derivative 
of Eq. ([5*5]) . we obtain d x b = (a — 3ga 2 )d x a. Using that d x b ~ a/A to eliminate 6, then, from 



a 



d x a = — — ^ (69) 

A(er — 3ga z ) 

we obtain that the discontinuity occurs when a 2 = a/3g. At that point b = era — ga 3 = 
-2a 3/2 /(3^3g), if we take the minus sign for a. As b is continuous through the discontinuity 
of APD (Fig. [6]), the value of a after the discontinuity can be calculated from 

ga 3 - aa + b = 0, (70) 



which gives a = 2ya/3g. Then, the conduction block occurs at a value of the period 
implicitly given by 

1 



DI min = r-A c - 2^/a/3g, with a = -f"{r - r c ), (71) 
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and at a point in space given by the position of the first singularity, or 

1 



x = -A(l-lnv / 3)r(r-r c ), (72) 



obtained solving Eq. ( 1691) with the initial condition a(0) = J a/ g. 

One should be very careful when considering this latter result, since the addition of 
diffusive effects changes the position of the node, and can make it travel. The reason is 
that, without diffusion, arbitrarily steep gradients can develop, that prevent the node from 
moving. Once diffusive effects are included, this is no longer the case, and the node travels, 
unless its motion is balanced by the drift term. In fact, there is a close analogy between the 
amplitude equation fl5"4"l) and the real Ginzburg-Landau equation that has been extensively 
studied in the context of phase transitions and front propagation j^sj. The dynamics is richer 
here because the integral term originating from dispersion causes a non-local interaction of 
the fronts separating two out-of-phase oscillating regions with the pacing end of the cable. 

This is easier to see assuming that both dispersion and the drift terms are small. Then, 
to first order, one recovers the Ginzburg-Landau equation 

rd t a = aa + ^ 2 d x a — ga 3 , (73) 



that has a stationary solution in the form of a front ao(x) = \J a / gtanh[-^/ a/(2£ 2 )(x — Xq)], 
connecting at x = xq regions oscillating out of phase. Then, to calculate the correction due 
to dispersion and asymmetric coupling, we can set a(x, t) = ao(x, xo(t)) + a\(x, t), which, 
assuming a\ small, results to first order in 

1 f x 

(cr + £ 2 dl — 3gal)a\ = rd t ao + wd x a + — ao(x')dx'. (74) 

A jo 

In order for the system to have a solution, the right hand side must be orthogonal to the 
left eigenvector of the linear operator, which in this case is simply d x ao since the operator 
is self-adjoint. Then, we obtain 

~ T ~1T I [SraoO^W+w 7 / [d x a Q (x')] 2 dx'+\ [ d x a (x r 
dt Jo Jo A Jo 

Evaluating the integrals, it is found that the motion of the node is given by 

dx 3 jl? 

T ^r = w -Ai^ x °- (76) 
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For arbitrary dispersion the expression for the motion becomes more complicated, but the 
effect is essentially the same. The gradient term corresponding to w makes the node move 
with constant velocity away from the pacing point, while the integral term provides a ramp 
that makes one phase of oscillation preferred over the other. When dispersion is weak, there 
is a point at which the two effects balance each other and the node relaxes to the equilibrium 
position 

wA [2a 

x ° 3 VF' 

which differs markedly from the prediction that neglects diffusive coupling given by Eq. 



E. Two-dimensional paced tissue 

Let us now generalize our formalism to the case of a square piece of tissue paced at one 
corner. The equivalent of Eq. ( 1431) is simply 

T d t a = aa-ga 3 -b-w(h- V)a + £ 2 V 2 a, (78) 

where n is a unit vector normal to the direction of propagation of the wavefront, and we 
impose non-flux boundary conditions on the boundaries of the tissue. 

Now the determination of the period at a given point would involve the integral along the 
path traveled by the wavefront. In particular, we obtain the condition for the oscillations in 
period: 

(n.V)6=~ (79) 

In general this is a complicated problem, since we have to know the normal to the wavefront 
at any given point. A major simplification occurs if we assume that the oscillations in APD 
do not affect the form of the propagating wavefront, which is the same as to say that the 
azimuthal variations of APD are small compared with the radial ones. Then, as long as the 
thickness of the wavefront is small compared with the size of the tissue, it can be assumed 
that, except for narrow boundary layers at the borders of the tissue, the propagation of the 
wavefront is perfectly circular. In this case, Eq. ( 1751) can be rewritten as: 

rr 

rd t a = aa — ga 3 — / adr' — w d r a + £ 2 V 2 a. (80) 
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In the far field limit, these equations reduce to Eq. (I54I) . with the node now becoming a 

circular line. We therefore expect to be a minimum tissue size for the node to form (Fie. 

fl 

[T3"a). as in the one- dimensional case. This was already observed in experiments [3j and 
numerical simulations of ionic models [7], where after a few beats, a nodal line was formed 
at the corner opposite the pacing point, and then moved towards it. Close to the pacing 
point, however, the motion of the node becomes nontrivial, because of curvature effects. 
Besides the term —wd r a that tends to make the nodal line move in the radial direction 
away from the pacing point, there is a contribution to the drift in the direction normal to 
the interface coming from the Laplacian, this being positive or negative depending on the 
curvature of the interface. As the nodal line is always perpendicular to the boundaries, its 
curvature depends on its position on the square. When the tissue is large enough so in the 
one-dimensional case the node would form in the region with positive curvature (close to 
the pacing point), this extra effect may be enough to make it travel (Fig. [T3|) . 



V. DISCUSSION 



T-wave alternans is thought to play a key role in the transition from normal heart rhythm 
to reentrant tachycardia and fibrillation. Alternans has been shown to induce spiral break- 



up 



361 ] . leading to a disordered spatio-temporal state that is usually identified with 



fibrillation. For this reason, it has been hypothesized that a steep restitution curve may be 
a possible determinant o: 



support this hypothesis 



the transition from VT to VF 



36j . Although some experiments 



371 ]. there is so far no conclusive evidence 38]. From a theoretical 



point of view, and despite earlier progress [8j, we still do not have a complete understanding 
on how alternans affects the stability of spirals. Numerical simulations show that break- 
up occurs when the amplitude of alternans oscillation grows enough to induce conduction 
blocks, which typically happens away from the spiral core. However, there is presently 
no theory that predicts the spatial distribution of alternans in a spiral wave. Topological 
arguments impose the existence of a nodal line (also termed line defect) 39|, with a jump in 
2tt in the phase of the oscillation, stretching from the core of the spiral to the boundaries. 
Thus, discordant alternans is always present, but what determines the motion of the nodal 
line, and whether it is important for spiral break-up is not entirely clear, except far from 
the spiral core where the dynamics is essentially one-dimensional. 
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In the present paper, we have considered the simpler cases of one and two-dimensional 
paced tissue and a circulating pulse in a ring geometry, and derived an equation for the 
spatio-temporal dynamics of small amplitude alternans in these states. Although conduction 
blocks fall out of the scope of the present theory, since it only considers small variations in 
conduction velocity, the results do shed light on where conduction blocks will occur. For 
that one can look at the distribution of APD and locate the places where it is larger than 
a critical value. To study the evolution after a conduction block has occurred, one has to 
resort again to simulations of the original equations. 

A circulating pulse in a ring has long been considered as a simplified model of anatomical 



reentry 



40j. For rings smaller than a critical value, oscillations in APD appear lCj. An 



understanding of these oscillations can be helpful in terminating anatomical reentry circuits. 
In this case discordant alternans is always present, and its wavelength is determined by the 
length of the ring, with the fastest growing mode having A ~ 2L. 

In the paced case, it has been shown that the gradients of repolarization created 
during discordant alternans offer a substrate for conduction block, leading to reentry and 

Tinna 

ventricular fibrillation. We have shown here that discordant wave patterns [3|, |5|, |6|, |7j result 
from a finite wavelength linear instability. Hence, their formation requires a minimum tissue 
size L min ~ A/4, required for at least one node to form. The value of L min that we measure 
in simulations of reaction-diffusion models are actually close to A/4 with A predicted by Eq. 
( l67l) and Eq. (|65l) . respectively (Table [IB . This lengthscale is similar in two-dimensional 
paced tissue, although in this case, curvature of the nodal line may affect the motion of this 
line. In both paced one- and two-dimensional tissue and the ring geometry, the onset of 
alternans in tissue is different than in a paced isolated cell because alternans, in tissue, is 
manifested as a wave and diffusive coupling tends to smooth out spatial gradients of APD. 

The paced and reentrant geometries studied in the present paper also form the basis to 
develop a theory of the dynamics of alternans in the presence of spiral waves, which has been 
the subject of both experimental and numerical studies 39]. Far from the core, the spiral is 
similar to the two-dimensional paced case. In the case where nodes move towards the pacing 
site in a one- dimensional cable, the nodal line should form a spiral with an opposite chirality 
as the propagating spiral wave front. This opposite chirality is imposed by the fact that the 
nodal line moves inward towards the core under the effect of a sufficiently steep, positively 
sloped, CV-restitution curve, and under the assumption that cellular alternans is driven by 
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APD restitution. The wavelength of the nodal line spiral far from the spiral core should 
be equal to the wavelength of discordant alternans in the one-dimensional paced case. This 
should match to the nodal line expected from one-dimensional reentry close to the spiral 
core. Then, the nodal line should extend straight out of the spiral core in the simplest limit 
of a constant wave speed where the wavelength of discordant alternans diverges far from the 
core. A theory for the spiral must, therefore, reduce to the two cases studied in this paper 
in the appropriate limits. The development of a theory of alternans in this case remains a 
fascinating task for future work. 

VI. CONCLUSIONS 

We have derived an equation that describes the spatio-temporal dynamics of small oscil- 
lation alternans in cardiac tissue. Our formulation is based on the restitution properties of 
the system, and takes also into account intercellular coupling, which is crucial in order to 
derive the correct threshold pacing rate and lengthscale of discordant alternans. For a sim- 
plified two-variable ionic model, we have been able to calculate these coefficients, and show 
that our reduced description agrees well with numerical simulations of the model. We have 
also considered the more realistic Noble model, and measured these coefficients numerically, 
stressing the generality of our formulation. We have applied our formulation to two different 
cases, a paced tissue and a circulating pulse in a ring. 

In the ring, the amplitude equation predicts that both dispersion and intercellular cou- 
pling affect the frequency of the motion, and the oscillations of APD are typically quasiperi- 
odic. There is a particular value of the length of the ring where these two effects balance 
each other, and the oscillations become periodic. For the paced case, the amplitude equa- 
tion predicts that discordant wave patterns result from a finite wavelength linear instability. 
Hence, their formation requires a minimum tissue size L m i n ~ A/4, required for at least one 
node to form. This lengthscale is intrinsic, in opposition to the circulating pulse, where it 
depends on the length of the ring. From the linear problem we deduce the possibility of two 
different patterns, depending on the parameters of the system, a standing wave pattern, or a 
traveling one. In the latter state, the amplitude of the oscillation grows exponentially away 
from the pacing point in a linear regime, which favors the induction of conduction blocks. 

Our formulation is based on the assumption that the simple map relationship (JTJ is 
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satisfied. The addition of memory effects would modify the coefficients of the amplitude 
equation, but not its general form, that is generic close to the period doubling bifurcation 
point. This modification, however, has important consequences for relating the genesis of 
discordant alternans to the underlying cell physiology. As will be shown elsewhere, one 
interesting effect of memory is that it makes the coefficient of the non-local term in Eq. 

depend both on the slope of the CV-restitution curve and the beat-to-beat dynamics. 
This opens the possibility to prevent the formation of discordant alternans by modifying 
the dynamics at a single-cell level to make this coefficient negative, as would be the case 
for a negatively sloped "supranormal" CV-restitution curve when the single cell dynamics 
is governed simply by APD-restitution. 

Even though the present theory does not take into account several effects present in 
real tissue (such as the three dimensional fiber geometry of an anatomical heart, the corre- 
sponding anisotropy in conduction velocity, spatial gradients of restitution properties, etc) 
it sheds light on the basic mechanisms and length scales that control the formation and evo- 
lution of discordant alternans. Moreover, it identifies a small subset of relevant parameters 
that control the formation of these arrhythmogenic patterns in experiments and numerical 
simulations of more complex ionic models. 
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APPENDIX A: ACTION POTENTIAL DURATION AND CONDUCTION VE- 
LOCITY OF THE TWO- VARIABLE MODEL 

In this appendix we show how to obtain the APD and CV-restitution curves for a prop- 
agating pulse in the two variable model. We will consider the limit e — » in which the 
sigmoidal becomes a step function S(V — V c ) — > Q(V — V c ). For a pulse propagating with 
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velocity c in the right direction, we can write dt = —cd x , and Eqs. 

DdlV + cd x V - 1/tq + h/r a = 



become: 



—cd x h = —h/r + 

Dd 2 x V + cd x V-V/(V c r )=0 
—cd x h = (1 — h)/r- 



when V > V r . 



when V <V r . 



(Al) 

(A2) 
(A3) 



and we will assume that, at x = we have the position of the wavefront given by V(0) = V c . 
The equations for the gate h can be integrated to give: 



h(x) = h e x/{cT+ \ when V > V c , 
h(x) = 1 - h ie x/icr ~\ when V < V c , 



(A4) 
(A5) 



where ho and hi are constants of integration. If we assume that APD >> r + (APD ~ 200 
ms, and r + = 12 ms in our simulations), then we can consider that h has decayed to zero by 
the end of the action potential, i.e., at the waveback of the previous pulse. Then h(cDT) = 0, 
and using Eq. flA5j) we obtain hi = exp (— DI/r_), and ho = h(0) = 1 — exp (— DI/r_). Next, 
we can solve the equation for the potential that, when V < V c becomes 



V{x) = Ae XlX + Be X2X , 



with 



Al,: 



1 



2D 



(A6) 



(A7) 



-ct^/c 2 + AD/(t V c ) . 

A solution in the wavefront (x > 0) must satisfy that V — ► when x is large and positive. 
Therefore it will correspond to a decaying mode: 



V(x > 0) = V c e Xlx = V c exp 



x 



_(_ c _^ + 4D /( To y e )) 



(A8) 



In the waveback (x < — APDc), on the other hand, we must impose that V — > when x is 
large and negative, corresponding to a growing mode in space. Then: 

'x + APDc 



V(x< -APDc) = v c e x ' 2{x+APDc) = V c 



exp 



2D 



-c+Jc* + 4D/(toV c )) 



(A9) 



The value of the potential for —APDc < x < will be given by the solution of 

Dd 2 x V + cd x V - 1/ro + h e x ^/r a = 0, 



(A10) 
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V C = A + B- - ; ± (A12) 



that can be written as: 

Wx) = Ae" ca;/I) + B + — - c ' r + ^V/cr +; for _ APDc < x < o, (All) 

cr D + c 2 t + r a 

where we recall that ho — 1 — exp (— DI/r_). Then, if we obtain the constants A and B 
imposing continuity of the potential and its derivative at i = 0, the same conditions at 
x = —APDc will give us the value of the APD and conduction velocity c. Imposing the bc's 
at x = and x = — cAPD we obtain the set of equations: 

c 2 t+ hp 

D + C 2 T + T a 

V c = -— + Ae c2APD / D + B- C M ^ e -APD/, + ; (A13) 

T D + C 2 T + T a 

V C X 2 = -^-Ae c2APD / D + — - " T + ^ e - APD ^. (A15) 

D CT D + C 2 T+ T a 

We have four equations for the four unknowns A, B, APD and c. These equations can be 
simplified provided that APD ^> r + and D jc <C cAPD. Then, we will define a new constant 
C = Ae c2AFD / D , and take the approximations exp(— c 2 APD/D) ~ 0, exp(— APD/r + ) ~ 0, 
from which we obtain 

APD 

V c = — + C+B, (A17) 

VA, = ± - J^^, (A18) 
cr D + c 2 t + r a 

V C \ 2 = -^C+—. (A19) 

D CT 

Eq. flAT8|) directly gives an implicit relation for the conduction velocity c in terms of the 
diastolic interval DI: 



ct_|_ 1 - e- m / T - cV c ( L AD \ 1 , . . 

£> + c 2 r+ r a 2D V V c 2 r y c y cr 

The value of the APD can be obtained from Eq. (1A17I) : 

APD = t (B + C -V c ), (A21) 

where 

r 2 T 2 h 

B = V c+J £p±-!±. (A22) 

D + C 2 T + T a 
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To obtain the constant C instead of directly using Eq. (1A19I) it is more convenient to take 
the difference of Eqs. f ]A18j) and (IA19I) from which we obtain 



C D + C Z T + T a 



Then, the APD becomes: 

AD 



APD = t 



^(l-e-^-)-Vjl 



(A24) 



T a V C 2 T V C 

To calculate the APD one has to solve then first Eq. (1A20I) for the conduction velocity. 
This can be avoided neglecting the term 4D / (c 2 TqV c ) in the square root, that introduces an 
error of 1 or 2 ms in the determination of the APD. This is equivalent to neglecting the 
effect of the outward repolarization current 1/tq during depolarization. In this way, CV and 
APD-restitution become decoupled, and we can write the final expression 

APD = Z±l£(i _ e - m/T -) - V C T . (A25) 
APPENDIX B: DERIVATION OF THE KERNEL 

Let us show how to calculate the kernel in Eq. (jHJ), and the coefficients w and £ 2 , for 
the two-variable model. This will allow us to verify the theory, as we can compare our 
analytical results with the simulations, and also gain some physical intuition on the origin of 
these coefficients. To derive the kernel, we start from the cable equation ([3]) (with I ext = 0). 
Our main goal is to quantify the effect of a gradient of APD on the restitution curve. For 
that, we will consider a single pulse propagating in this gradient with speed c. Interpreting 
the cable equation as a diffusion equation where I ion is a spatially distributed source, Eq. (j3J) 
can be formally inverted, expressing the transmembrane potential V in terms of the Green 
function of the diffusion operator: 

r (x—x') 2 \ 

Vix t )=— c dt' r dx' —^^t± )±i. u> t ') (Bi) 

V[X ^ C m J-oo loo [47rD(t-t')] 1/2 
where we have assumed translational invariance, as is the case in the ring. 

Now, if we are measuring the action potential duration at a given value of the transmem- 
brane potential V c , then, by definition 

V c = V{x,x/c + A(x)), (B2) 
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where x/c is the time it took the pulse to arrive at position x and, for simplicity we introduce 
the notation A(x) = APD (a;). Substituting this into Eq. fIBll) . we obtain 

V -— f X/C+A(X) dt> r dx> eXp {~^(^X)-^)) j m fB3 ) 

This defines an implicit integral equation for the action potential duration. 

In the general case it is not clear how to invert this equation. We can do it, however, for 
the simpler case of the two variable model, where the action potential adopts a triangular 
form. Integrating Eq. (jSJ) for the gate variable h (in the limit e — > 0), for V > V c , we obtain 
h = hoe~ t / T+ , where ho = 1 — e~ m / T - can be obtained integrating Eq. ((Sj) for V < V c , and 
assuming r + APD, so by the end of the previous action potential the gate is completely 
closed (h = 0). 

Then, when V > V c , and in the absence of propagation, the equation for the temporal 
evolution of the transmembrane potential becomes: 

V = -{l-e~ m / T -)e- t/T + --. (B4) 

Ta TO 

When t + ro we can assume that depolarization occurs instantaneously, so we can write 
the current in the form 

kon(t)/C m = -6(t)I + 9(t - APD)/ + 5(t) J[DI], (B5) 

where 9(t) is the standard Heaviside step function, and 5(t) the Dirac delta function. The 
former equation means that, after an excitation at t — 0, the voltage takes the maximum 
value V max = J[Dl] and then decreases linearly in time V = V max —It, with the identifications 
J[DI] = r+(l - e- DI /-)/ Ta and / = l/r . 

Integrating Eq. (1B4I) . under the approximation r + ^ tq we obtain 

V{t) = ^±(1 - e- m ' T -) -- = V max - -. (B6) 

At t = APD, V(APD) = Vc, and the value of the APD becomes APD = (V max - V c )r , or 

APD n+i = (j|pi»] _ v c )// = Z±Z£(i _ e" DI / T -) - Kro = /(DP), (B7) 

which gives us the APD-restitution. This is the same expression that we obtained in the 
previous appendix (cf. Eq. ( 1A25I) ). neglecting the term 4D / (c 2 ToV c ) . 
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For a propagating pulse in a gradient of DI, we can again use expression (1B5I) . but taking 
into account that the excitation now occurs when the pulse arrives (i.e. at time t = x/c), 
and that the diastolic interval depends on space. Then 



I ion (x, t)/C m = -6(t - x/c)I + 6(t - x/c - A(x))I + 8(t - x/c)J[Dl(x)}. 



(B8) 



It is important to emphasize that V max = J[DI] is different for an isolated cell than for a 
cable because of diffusive coupling. Therefore J[DI] in Eq. ( 1B8I) denotes the peak value of 
the voltage for a propagated action potential, as opposed to a stimulated cell in Eq. (1B5I) . 
In contrast, the local repolarizing current I is the same in both cases. 

Substituting this expression into (1B3|) we can now split the integral into three parts 
V c = Vi + Vn + Vin. The first part is 



V, 



x/c+A(x) roc exp I 

dt' dx' [ 



(x-x') 2 



4D(x/c+A(x)- 



r x/c+A(x) ret' exp I 

I dt' dx' -i 

J-oo J-oo \4-KD(a 



[4ttD{x/c + A{x) -t')} 1 / 

(x-x') 2 \ 

t')S 



4D(x/c+A(x)-t 



This gives 



rx/c+A(x) I 

Vj = I J ^ dt'- {1 + Erf 



[AttD{x/c + A(x) -t')} 1 / 2 ' 
ct' — X 



4D(x/c + A(x)-t') 



(B9) 



(BIO) 



Making the approximation 

1 



1 + Erf 



4D(A(x) - y/c) 



(Bll) 



valid when D/c < (DA) 1/2 < Ac, we obtain 

fx/c+A(x) 



Vi 



dt'9(ct' -x) = IA(x). 



(B12) 



As long as the length scales associated with intercellular coupling are small compared with 
the distance the pulse travels during an APD (so the coupling with other cells is negligible), 
this integral gives us the change in voltage during repolarization at a given point, during 
the time of one APD. 

The second part becomes: 



V, 



11 



x/c+A(x) poo exp 

dt' / dx' 



(x-x 1 ) 2 



4D(x/c+A(x) 



[AttD(x/c + A(x 



^A_9(t'-x'/c-A(x'))I 



(B13) 
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To be able to solve for Vji we make the approximation A[x) ~ A c . Then: 

/ „ , , » N f (X-X') 2 1 

p/c+A c rc(t'-A c ) ex P 1 ~ 4D(Wc+A f 



/ / dt'- ll + Erf 



c 



:Jx/c+A c - t' 



(B14) 



\/AD 

Making the change of variable y = c 2 (x/c + A c — £') / (4D) we obtain: 

2D r°° D 
V H ~—lj Q dy[l + Erf(-^)} = -I (B15) 

Clearly this term is much smaller that the previous one (Vj-j/Vj ~ (D/c)/(cA c ) <C 1) and 
we will neglect it. 
The last term is 

Vttt - r' C+A(X) dt' r dx' ^("^K)) J\m(x')]8(t'-x'/c) 
VlH - J-oo d J-oo dX [4ttD(x/c + A(x) -t')] 1/2 [ 1 )l [ ' ] 

X/C+A{X) dt' ^{"^ ( ^ k^h cJ [DI(^)]. (B16) 
[4nD(x/c + A(x) -t')} 1 ' 2 [ v n K ' 

Making the change of variable t' = (x + y)j c, we finally obtain 

VlII = L dv [^D{A{x)-y/c)]V* J[m{x + v)] - (B17) 
And this integral gives us the correction to the maximum value of the APD coming from 
the intercellular coupling, as Vm ~ V max . 

Adding Vi and Vm, and solving for A(x) in Eq. ( 1B12I) . we obtain 

> Z [ w^k m{x+v)l (B18) 

where A(x) is the APD following the diastolic interval DI, so A(x) = APD n+1 and DI = 
DI n , and we have extended the limits of the integral to infinity, assuming a rapid decay of 
the kernel. Note that A(x) is still in the kernel, so the former equation gives an implicit 
expression for the APD. As we are interested in the regime close to the onset of period 
doubling, we can approximate the APD within the integral as A(x) ~ A c . The deviations 
from the critical value would give nonlinear gradient terms in the oscillations of APD, that 
we assume to be of higher order. Then, our final expression is: 

APD- (*) = £ dv ^t0E^ m r (x + ,)]. (B19) 
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Once we know the expression for the kernel, the coefficients w and £ 2 are easy to obtain. It is 
useful to rewrite the kernel using the change of variable y = 2(D Ac) 1 / 2 'z. Then, expanding 
the kernel to first order in the coefficient 2(D A C Y^ 2 / (A c c), again assumed to be small, we 
get 

1 r°° o 1 FT) 

APD n+1 (x) = — / dze' z [1 + -J — (z-2z 3 )}f[DI n (x + 2(DA c ) 1 / 2 z)}. (B20) 

\ 7T J-oo C V A c 



We can expand the restitution curve around its value at a given point 

f[DF(x + y)] = f [DI» + yd x Dl n (x) + \y 2 d 2 x DF(x) + ■■■ 

= f[DI n (x)} + 2z{DA c ) 1 ' 2 fd x m n {x) + 2z 2 DA c fd 2 x m n {x) 

+2z 2 DA c f"[d x m n (x)} 2 + --- (B21) 

Introducing this expansion in Eq. (1B20I) . the asymmetrical part results 
2D 1 r°° i 2D 

±—faj>F(x)-= / dze- z \z 2 - 2z A ) = fd x Dr(x), (B22) 

C V7T J-oo C 

and the symmetrical one 

2DA C {fd 2 x Dr(x) + f"[d x DI n (x)] 2 } dze~ z2 z 2 = DA C {f'd 2 x DF(x) + f"[d x Dr(x)] 2 } 

(B23) 

In this last expresion we will neglect the term [<9 x DI n (a;)] 2 since, close to the onset of alter- 
nans, it is higher order with respect to the term d 2 Dl n (x) (their ratio is ~ (/' / f")5D, where 
5D is the amplitude of the oscillations in DI(x)). However, if the system is far from onset, 
it maybe necessary to take it into account, if we want to obtain quantitative results. 
Then, we have 

2D 

APD n+1 (x) = /[DI n (x)l fd x Dl n (x) + DA c fd 2 m n {x) + ■■■ (B24) 

c 

This expression is the same as Eq. (HOI) . Identifying terms we obtain the values w = 2D/c, 
£ 2 = D APD C . For more complicated ionic models these coefficients have to be calculated 
numerically following the procedure outlined in Sec. IV. A. 
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FIG. 1: a) Trajectories in the phase plane (V, h) for the two-variable model [Eqs. (HJ)-©], obtained 
by simulating V = —Iion/C m with an activation interval r = 400 ms and the parameters given in 
the text. The dashed lines denote the nullclines h = and V = 0. b) Time evolution of the voltage 
V (solid line) and the gate variable h (dashed line) during such trajectory. 
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FIG. 2: APD- and CV-restitution curves corresponding to the Noble [a) and c)] and two-variable 
model [b) and d)]. The solid lines correspond to the restitution curves obtained numerically using 
a S1-S2 protocol. In b) and d) the dotted lines correspond to the theoretical curves given by 
Eqs. (|A25j) and (|A20p . The dotted-dashed line in d) corresponds to simulations with a finer grid 
dx=0.001, and dt=0.001, which agree very well with the theoretical predictions, and indicate that, 
for the usual parameters in the simulations, there is an error of about a 4% in the determination 
of the maximum conduction velocity. 
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FIG. 3: Left panel: distribution of APD at consecutive beats, obtained solving Eq. §3§ in a ring 
of length L = 5 cm, for the two-variable ionic model. Right panel: evolution in time of the value 
of the APD at x = L/2 (top), and position of the node (bottom), defined as the point xq at which 
APD(xq) = (APD max + APD min )/2. 
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FIG. 4: Bifurcation diagram obtained from simulations of the two- variable model in the ring. The 
bifurcation is supercritical, with a critical ring length of L c ~ 5.11 cm. 



47 




5 10 15 20 

L (cm) 




5 10 15 

L (cm) 



FIG. 5: Distribution of the action potential duration in a long strand of tissue obtained simulating 
Eq. ([3]) for the Noble a), b), and c) and two- variable models, d), with r = 258ms, and r = 290ms, 
respectively, and several tissue lengths. The solid and dashed lines correspond to two consecutive 
beats, while the dotted line in d) represents the APD ten beats later. Thus, in d) the pattern is 
traveling towards the pacing point and the amplitude of alternans grows away from the pacing site 
as predicted theoretically. 
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FIG. 6: Upper panel: Oscillations in action potential duration and period, obtained simulating 
Eqs. (|53p and (|55p . with the parameters of the two- variable model, and r = 290ms. Starting with 
a constant value of a, we show in dashed lines the evolution of the system at times t = 0, 10r, 50r, 
and 500t. The final state is shown in solid lines. Below we sketch the bifurcation diagram at 
several points along the cable for this final state. The dashed line denotes the value of T(0) = r 
and the filled circle the value of the amplitude of oscillation in APD in each case. As we go along 
the tissue b becomes positive and the pitchfork (dotted line) becomes and imperfect bifurcation 
[cf. Eq. (|55p ] (a). At point (b) the saddle-node is at T = T(0). There is a jump to the other phase 
and b starts to decrease (c). At (d) 6 = and a perfect pitchfork bifurcation is recovered, after 
which the branches switch place (e). Finally, the saddle-node reaches T = T(0) and the process 
starts again. 
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DI (ms) 

FIG. 7: a) Distribution of the action potential potential at two consecutive beats in a long strand 
of tissue obtained simulating Eq. ([3]) for the Noble model, with r = 258ms. b) Restitution 
curves for the Noble model. The dotted line corresponds to that obtained using a S1-S2 protocol, 
while the solid and dashed lines are the restitution curves corresponding to the two beats in 
a). The difference between these two latter curves at a point x* corresponding to a node in DI, 
AAPDi = APD n+1 (x*) — APD n (x*), gives the change in APD corresponding to a negative or 
positive slope. The value AAPD2 is obtained as the difference between the maximum value of 
the APD in a), and the value of APD for the same DI obtained using a S1-S2 protocol, where no 
gradients are present. 
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FIG. 8: Comparison of the results obtained simulating the Noble model (solid line) and the am- 
plitude equation (|5H) (dashed line) , with the coefficients given in Table HH 




FIG. 9: Convective vs absolute instability. When the instability is convective (as in a)), an initial 
localized perturbation is advected as it grows (here £3 > ti > ii), so at a given point in space it 
decays. The instability becomes absolute when the wavepacket grows at any given point in space 
(b). This is signaled by a vanishing group velocity. It should be noted that the phase velocity, 
however, does not vanish and, in general, it does not have to be in the same direction as the group 
velocity. 
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FIG. 10: (Color online). Top panels: comparison between the spectra obtained solving the linear 
eigenvalue problem given by Eq. (|66p (circles) and the analytical predictions (diamonds and stars), 
for the parameters of Noble model (left) at a pacing period of r = 258 ms and L = 20 cm, and 
the two- variable model (right), with r = 290 ms and L = 40 cm. Bottom panels: dependence of 
growth rate and frequency on the real wavenumber. For the linear eigenvalue problem (circles), 
the wavenumber has been calculated from the number of nodes n of the eigenmodes k r ~ nir/L. 
The star corresponds to the stationary mode, with k r = (u>A) -1 / 2 , Qi = 0, and f2 r = a — £ 2 /(wA). 
The diamonds correspond to the results for the standing waves, obtained solving Eq. (|62p . For 
the parameters of Noble (left) the first mode to bifurcate is an stationary mode, corresponding to 
the isolated eigenvalue. For the two- variable model (right), there is a continuous branch of modes 
that bifurcates first. 
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FIG. 11: Space-time plots of a obtained by simulations of Eq. (|54p for parameters of the two- 
variable model, showing absolutely unstable (left and r = 295 ms) and convectively unstable (right 
and r = 298 ms ) wave patterns. The crosses denote the positions of the nodes (a = 0). 
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FIG. 12: Stability diagram of the Noble (left panel) and two-variable (right panel) cable models 
with domains of no-alternans (open circles), concordant alternans (filled squares), and discordant 
alternans (filled diamonds); conduction blocks form at smaller r not shown here. Boundaries 
between the same domains obtained by simulations of the amplitude equation ([5^1 ) are shown by 
solid lines. The dashed line denotes the bifurcation period for alternans predicted by the map of 
Eq. ©■ 
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FIG. 13: Simulations of Eq. (|80p for the parameters of the Noble model, with r = 258 ms and a) 
L = 2.5 cm, b) L = 3.5 cm. The solid lines represent the position of the node. The lines in b) are 
drawn every ten beats. As the size of the tissue increases, a node forms at the corner opposite the 
pacing point (a), which above a certain tissue size begins to travel, due to curvature effects (b). 
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